#include <cstdio>
#include <string>
#include <iomanip>
#include <algorithm>
#include <iostream>
#include <istream>
#include <stdexcept>
#include <algorithm>
#include <vector>
#include <boost/algorithm/string/split.hpp>
#include "is_equal_to.h"
#include "ultoa.h"
#include "perlutil.h"
#include "blast_results.h"
#include "bcpp_read_seq.h"
using std::cerr;
using std::string;
using std::istream;
using std::ostream;
using std::vector;


//#include "alignlib/alignlib.h"
//#include "alignlib/Alignandum.h"
//#include "alignlib/AlignataIterator.h"
//#include "alignlib/Alignata.h"
//#include "alignlib/AlignException.h"
//#include "alignlib/HelpersAlignata.h"
//#include "alignlib/Translator.h"
//#include "alignlib/HelpersTranslator.h"
//#include "alignlib/Alignatum.h"
//#include "alignlib/HelpersAlignatum.h"
//#include "alignlib/Alignator.h"
//#include "alignlib/HelpersAlignator.h"
//#include "alignlib/SubstitutionMatrix.h"
//#include "alignlib/HelpersSubstitutionMatrix.h"





/* Used by the callback function. */

namespace
{
	char* initialize_table(char* table)
	{
		const string aminoacids = " ABCDEFGHIKLMNPQRSTVWXY";
		for (char i = 0, len = (char)aminoacids.length(); i < len; ++i)
		{
			table[aminoacids[i]] = i;
		}
		return table;
	}

	char* initialize_reverse_table(char* table)
	{
		const string aminoacids = " ABCDEFGHIKLMNPQRSTVWXY";
		for (char i = 0, len = (char)aminoacids.length(); i < len; ++i)
		{
			table[i] = aminoacids[i];
		}
		return table;
	}
}

namespace bcpp_align
{
	void bl_prepare_sequence_str(string& str)
	{
		static char convert_table[256];
		static char* p_convert_buf = initialize_table(convert_table);
		string::iterator i = str.begin();string::iterator end = str.end();
		for (; i != end; ++i)
			*i = p_convert_buf[*i];
	}

	string bl_retrieve_string_from_prepared_sequence(const string& param)
	{
		static char convert_table[256];
		static char* p_convert_buf = initialize_reverse_table(convert_table);
		string str = param;
		string::iterator i = str.begin();string::iterator end = str.end();
		for (; i != end; ++i)
			*i = p_convert_buf[*i];
		return str;
	}


	void cdna_gaps_to_prot_gaps(std::vector<t_gap>& cdna_gaps)
	{
		for (unsigned i = 0; i < cdna_gaps.size(); ++i)
		{
			cdna_gaps[i].len /= 3;
			cdna_gaps[i].beg /= 3;
		}
	}


	void get_gaps_from_cdna_sequence(const string& sequence, vector<t_gap>& gaps, char gap_char)
	{
		get_gaps_from_prot_sequence(sequence, gaps, gap_char);
		cdna_gaps_to_prot_gaps(gaps);
	}

	void get_gaps_from_prot_sequence(const string& sequence, vector<t_gap>& gaps, char gap_char)
	{
		gaps.clear();

		//save the number and position of the aligning characters for each sequence
		unsigned cummulative_offset = 0;
		string::const_iterator begin = sequence.begin();
		string::const_iterator pos = begin;
		string::const_iterator end = sequence.end();
		while (pos != end)
		{
			// find gap char
			pos = find(pos, end, gap_char);
			if (pos == end)
				break;

			string::const_iterator gap_start = pos;

			// find one past stretch of gaps
			while (++pos != end)
			{
				if (*pos != gap_char)
					break;
			}

			unsigned len = pos - gap_start;
			gaps.push_back(t_gap(gap_start - begin - cummulative_offset, len));
			cummulative_offset += len;
		}
	}


	inline void gap_to_str(const t_gap& gap, string& str)
	{
		ultoa(gap.beg, str) += ':';
		ultoa(gap.len, str);
	}


	string gap_array_to_str(const vector<t_gap>& gaps)
	{
		string gap_str;
		if (!gaps.size())
			return string();
		gap_to_str(gaps[0], gap_str);

		for (unsigned i = 1; i < gaps.size(); ++i)
		{
			gap_str += '|';
			gap_to_str(gaps[i], gap_str);
		}
		return gap_str;

	}


	bool gap_str_to_array(const string& gaps_str, vector<t_gap>& gaps)
	{
		// empty string encodes no gaps
		if (! gaps_str.length())
			return true;


		// split string of gaps
		vector<string> gap_segments;
		boost::algorithm::split( gap_segments, gaps_str, is_equal_to('|'), 
								 boost::token_compress_on);
		//perl_split(gaps_str,  "|",  perlutil::eCOUNT_EMPTY, gap_segments);
		for (unsigned i = 0; i < gap_segments.size(); ++i)
		{
			string::size_type pos = gap_segments[i].find(':');
			if (pos == string::npos)
				return false;
			// gaps should not come after a protein or be humungously long
			unsigned beg = atoi(gap_segments[i].c_str());
			unsigned len = atoi(gap_segments[i].c_str() + pos + 1);
			if (len > 80000)
				return false;
			gaps.push_back(t_gap(beg,  len));
		}
		return true;
	}

	bool expand_gaps_prot(const vector<t_gap>& gaps, string& prot)
	{
		for (int i = gaps.size()-1; i >= 0; --i)
		{
			if (gaps[i].beg > prot.length())
				return false;
			assert(gaps[i].len < 80000);
			prot.insert(gaps[i].beg, string(gaps[i].len,  '-'));
		}
		return true;
	}

	bool expand_gaps_cdna(const vector<t_gap>& gaps, string& cdna)
	{
		for (int i = gaps.size()-1; i >= 0; --i)
		{
			if (gaps[i].beg * 3 > cdna.length())
				return false;
			assert(gaps[i].len * 3  < 80000);
			cdna.insert(gaps[i].beg * 3, string(gaps[i].len * 3,  '-'));
		}
		return true;
	}

	ostream& operator << (ostream& ostrm, const t_align_results& results)
	{
		ostrm
		<< results.align_beg                            << "\t"
		<< results.align_len                            << "\t"
		<< results.total_len                            << "\t"
		<< std::setprecision(3)
		<< results.aligned_residues * 1.0 / 
		std::max(results.len1, results.len2)    << "\t"
		<< std::setprecision(5)
		<< results.e__value                             << "\t"
		<< results.len1                                 << "\t"
		<< results.len2                                 << "\t"
		<< results.aligned_residues                     << "\t"
		<< results.percent_id                           << "\t"
		<< results.bitscore                             << "\t"
		<< gap_array_to_str(results.gaps1)              << "\t"
		<< gap_array_to_str(results.gaps2);


		return ostrm;
	}


	char* initialise_valid_gaps()
	{
		static char valid_gaps[256];
		std::fill(valid_gaps,  valid_gaps + sizeof(valid_gaps),  '\x1');
		for (unsigned i = '0'; i <= '9'; ++i)
			valid_gaps[i] = 0;
		valid_gaps[static_cast<unsigned>('|')] = 0;
		valid_gaps[static_cast<unsigned>(':')] = 0;
		return valid_gaps;
	}
	string read_gap_str(istream& istrm)
	{
		// ignore first tab
		// leave trailing whitespace

		static char* ptable = initialise_valid_gaps();
		// ignore first tab
		char tabchar = istrm.get();
		if (tabchar != '\t')
			throw std::runtime_error("Alignment gaps string in invalid format.");
		string gaps_str;
		char c;
		while ((c = istrm.peek()) != EOF)
		{
			if (ptable[static_cast<unsigned>(c)])
				break;
			gaps_str += istrm.get();
		}
		return gaps_str;
	}
	istream& operator >> (istream& istrm, t_align_results& results)
	{
		double dummy;

		istrm
		>>  results.align_beg;
		if (!istrm)
			return istrm;

		istrm   >>  results.align_len
		>>  results.total_len
		// coverage not stored
		>>  dummy
		>>  results.e__value
		>>  results.len1
		>>  results.len2
		>>  results.aligned_residues
		>>  results.percent_id
		>>  results.bitscore;



		// get gaps from strings
		string gaps1 = read_gap_str(istrm);
		results.gaps1.clear();
		if (gaps1.length() && istrm)
			gap_str_to_array(gaps1, results.gaps1);

		string gaps2 = read_gap_str(istrm);
		results.gaps2.clear();
		if (gaps2.length() && istrm)
			gap_str_to_array(gaps2, results.gaps2);
		return istrm;
	}

	string t_align_results::get_sql_column_names()
	{
		return
		"align_start,"
		"align_len,"
		"aligned_positions,"
		"coverage,"
		"e_value,"
		"len1,"
		"len2,"
		"aligned_residues,"
		"percent_id,"
		"bitscore,"
		"gaps1,"
		"gaps2";
	}

	string t_align_results::column_names()
	{
		return
		"start\t"
		"len\t"
		"positions\t"
		"coverage\t"
		"e_value\t"
		"len1\t"
		"len2\t"
		"residues\t"
		"percent_id\t"
		"bitscore\t"
		"gaps1\t"
		"gaps2"
		;
	}

	string t_align_results::column_name_underscores()
	{
		return
		"-----\t"
		"---\t"
		"---------\t"
		"--------\t"
		"-------\t"
		"----\t"
		"----\t"
		"--------\t"
		"----------\t"
		"--------\t"
		"-----\t"
		"-----"
		;
	}

	string t_align_results::column_names2()
	{
		return
		"start1\t"
		"end1\t"
		"start2\t"
		"end2\t"
		"positions\t"
		"coverage\t"
		"e_value\t"
		"len1\t"
		"len2\t"
		"residues\t"
		"percent_id\t"
		"bitscore\t"
		"gaps1\t"
		"gaps2"
		;
	}

	string t_align_results::column_name_underscores2()
	{
		return
		"------\t"
		"----\t"
		"------\t"
		"----\t"
		"---------\t"
		"--------\t"
		"-_-----\t"
		"----\t"
		"----\t"
		"--------\t"
		"-------_--\t"
		"--------\t"
		"-----\t"
		"-----"
		;
	}

	ostream& operator << (ostream& ostrm, const t_debug_align_results& debug)
	{
		ostrm << "e-value           = " <<debug.blast_results.e__value                << "\n";
		ostrm << "bitscore          = " <<debug.blast_results.bitscore                << "\n";
		ostrm << "% identity        = " <<debug.blast_results.percent_id              << "\n";
		ostrm << "gaps1             = " <<gap_array_to_str(debug.blast_results.gaps1) << "\n";
		ostrm << "gaps2             = " <<gap_array_to_str(debug.blast_results.gaps2) << "\n";
		ostrm << "align_beg         = " <<debug.blast_results.align_beg               << "\n";
		ostrm << "align_len         = " <<debug.blast_results.align_len               << "\n";
		ostrm << "total_len         = " <<debug.blast_results.total_len               << "\n";
		ostrm << "length seq1       = " <<debug.blast_results.len1                    << "\n";
		ostrm << "length seq2       = " <<debug.blast_results.len2                    << "\n";
		ostrm << "aligned residues  = " <<debug.blast_results.aligned_residues        << "\n";
		ostrm << "error_message     = " <<debug.blast_results.error_message           << "\n";
		return ostrm;
	}


/*
ostream& output_no_alignment (ostream& ostrm)
{
	ostrm << "-1\t"             // bl2seq_bitscore
			 "\\N\t"            // bl2seq_align_start
			 "\\N\t"            // bl2seq_align_len
			 "\\N\t"            // bl2seq_total_len
			 "\\N\t"            // bl2seq_coverage
			 "-1\t"             // bl2seq_e_value
			 "\\N\t"            // aln_gaps1
			 "\\N\t"            // aln_gaps2
			 "\\N\t"            // aln_gaps2
			 "\\N\t"            // len1
			 "\\N\t"            // len2
			 "0";				// aligned residues
	return ostrm;

}
*/

/*
void bl_prepare_all_sequences(BioCPP::t_collect_sequences& seq_db)
{
	for (unsigned sz = seq_db.size(), i = 0; i < sz; ++i )
		bl_prepare_sequence_str(seq_db.sequence(i));
}

*/
	namespace
	{
		void count_cdna_aligned(const string& seq1, const string& seq2,
								unsigned& align_positions,
								unsigned& align_residues)
		{
//	string cdna_index1 = seq1;
//	string cdna_index2 = seq2;
//	validate_cdna_to_index(cdna_index1);
//	validate_cdna_to_index(cdna_index2);
//	unsigned cnt_identical = 0;

			string::const_iterator i1 = seq1.begin();
			string::const_iterator i2 = seq2.begin();
			string::const_iterator end = seq1.end();
			for (; i1 != end; i1 += 3, i2 +=3)
			{
				bool gap1= (*(i1 + 0) == '-' ||
							*(i1 + 1) == '-' ||
							*(i1 + 2) == '-');
				bool gap2= (*(i2 + 0) == '-' ||
							*(i2 + 1) == '-' ||
							*(i2 + 2) == '-');
				// gaps in both, ignore
				if (gap1 && gap2)
					continue;

				// else not gap in at least one, must be an aligned position
				++align_positions;

				// gaps in neither: aligned residue
				if (!gap1 && !gap2)
					++align_residues;

			}
		}
		void count_prot_aligned(const string& seq1, const string& seq2,
								unsigned& align_positions,
								unsigned& align_residues)
		{
			string::const_iterator i1 = seq1.begin();
			string::const_iterator i2 = seq2.begin();
			string::const_iterator end = seq1.end();
			for (; i1 != end; i1 ++, i2 ++)
			{
				bool gap1= *i1 == '-';
				bool gap2= *i2 == '-';
				// gaps in both, ignore
				if (gap1 && gap2)
					continue;

				// else not gap in at least one, must be an aligned position
				++align_positions;

				// gaps in neither: aligned residue
				if (!gap1 && !gap2)
					++align_residues;

			}
		}

		void get_align_results_from_sequence_pairs(t_align_results& align_result,
												   string seq1,
												   string seq2, bool cdna)
		{
			unsigned align_positions = 0;
			unsigned align_residues = 0;
			if (cdna)
				count_cdna_aligned(seq1, seq2, align_positions, align_residues);
			else
				count_prot_aligned(seq1, seq2, align_positions, align_residues);

			align_result.align_len          =
			align_result.total_len          = align_positions;
			align_result.aligned_residues   = align_residues;

			// get gaps
			if (cdna)
			{
				get_gaps_from_cdna_sequence(seq1, align_result.gaps1);
				get_gaps_from_cdna_sequence(seq2, align_result.gaps2);
			}
			else
			{
				get_gaps_from_prot_sequence(seq1, align_result.gaps1);
				get_gaps_from_prot_sequence(seq2, align_result.gaps2);
			}

			// calculate sequence lengths
			seq1.erase(std::remove(seq1.begin(), seq1.end(),  '-'), seq1.end());
			seq2.erase(std::remove(seq2.begin(), seq2.end(),  '-'), seq2.end());
			align_result.len1 = seq1.length();
			align_result.len2 = seq2.length();

			if (cdna)
			{
				align_result.len1 /= 3;
				align_result.len2 /= 3;
			}
		}
	}

	void get_align_results_from_cdna_sequence_pairs(t_align_results& align_result,
													const string& seq1,
													const string& seq2)
	{
		get_align_results_from_sequence_pairs(align_result, seq1, seq2, true);
	}

	void get_align_results_from_prot_sequence_pairs(t_align_results& align_result,
													const string& seq1,
													const string& seq2)
	{
		get_align_results_from_sequence_pairs(align_result, seq1, seq2, false);
	}



//________________________________________________________________________________________

//		insert_gaps_into_sequence_pair

//________________________________________________________________________________________
	void insert_gaps_into_sequence_pair(    const t_align_results&      pair_data,
											string&                     sequence1,
											string&                     sequence2,
											bool                        cdna,
											string&                     acc1,
											string&                     acc2)
	{
		unsigned codon = cdna ? 3 : 1;
		if (!pair_data.align_len)
			throw std::runtime_error(acc1 + " vs " + acc2 + ": Zero length alignment!");


		if (!sequence1.length())
			throw std::runtime_error("Sequence for " + acc1 + " is missing.");
		if (!sequence2.length())
			throw std::runtime_error("Sequence for " + acc2 + " is missing.");

		if (cdna)
		{
			if (!expand_gaps_cdna(pair_data.gaps1, sequence1) ||
				!expand_gaps_cdna(pair_data.gaps2, sequence2))
				throw std::runtime_error("Alignment gaps for " +
										 acc1 + " vs " + acc2 + 
										 " do not match sequence.");
		}
		else
		{
			if (!expand_gaps_prot(pair_data.gaps1, sequence1) ||
				!expand_gaps_prot(pair_data.gaps2, sequence2))
				throw std::runtime_error("Alignment gaps for " +
										 acc1 + " vs " + acc2 + 
										 " do not match sequence.");

		}
		if (sequence1.length() != sequence2.length())
			throw std::runtime_error("Sequence lengths for " +
									 acc1 + " and " +  acc2 + 
									 " are different after inserting gaps.");

		if (pair_data.align_beg * codon > sequence1.length())
			throw std::runtime_error("Alignment start position too large for " +
									 acc1 + " vs " +  acc2);

		// remove non-aligned region
		sequence1.erase(0, pair_data.align_beg * codon);
		sequence2.erase(0, pair_data.align_beg * codon);

		if (pair_data.align_len * codon > sequence1.length())
			throw std::runtime_error("Alignment length too large for " +
									 acc1 + " vs " +  acc2);

		sequence1.erase(pair_data.align_len * codon);
		sequence2.erase(pair_data.align_len * codon);
		return;

	}
	


//0123456789012345678901234567890123456789
//0123456789012345678901234567890123456789
//---AAA--AAAAAA--AAAAAAA
//   012  345678  90123456
//AAAAAAAAAAA--AAAAAAA---
//01234567890  1234567   
//    |<---------->|
//align_start=4
//align_len=14 (18 - 4)
//gaps1 = 0:3|3:2|9:2
//gaps2 = 11:2|18:3
//
//Andreas 
//start1 = 1
//end1 = 12
//start2 = 4
//end2 = 16
//
//gaps1 = +2-2+6-2+2
//gaps1 = +7-2+5
//
//
//gaps1 = 0:(start2 - start1)
//Running pos = start1		1
//align_len = gaps1+start1	4
//Running pos += +length; 	3
//align_len += (+length)		6   
//Gap length = (-length);		2
//align_len += (-length)		8   
//3:2                     	
//Running pos += +length; 	9
//align_len += (+length)		14  
//Gap length = (-length);		2
//align_len += (-length)		16  
//9:2
//Running pos += +length; 	11
//align_len += (+length)		18
//Unaligned = (len1 16 - Running_pos 11)	5
//Running pos += Unaligned; 	16
//align_len += Unaligned		23
//	
//gaps2 = 0:(start1 - start2)	0
//Running pos = start2		4
//align_len = gaps2+start2	4
//Running pos += +length; 	11
//align_len += (+length)		11   
//Gap length = (-length);		2
//align_len += (-length)		13   
//11:2                     	
//Running pos += +length; 	16
//align_len += (+length)		18  
//Unaligned = (len2 18 - Running_pos 16) 2
//Running pos += Unaligned; 	18
//align_len += Unaligned		20
//end gap len = other.align_len - align_len
//Running pos
	



//--------------------------------------------------------------------------------------------------
/** write an alignment in compressed form into two streams
	using namespace AlignLib;
	void alignata_to_blast_results( const AlignLib::Alignata& src, 
									const string& seq1,
									const string& seq2,
									t_align_results& res)
	{
		// sanity checks
		if (src.isEmpty())
			return;

		unsigned curpos1 = src.getRowFrom() - 1;
		unsigned curpos2 = src.getColFrom() - 1;
		res.len1 = seq1.length();
		res.len2 = seq2.length();
		res.align_beg = curpos1;

		// initial gaps
		if (curpos2 > curpos1)
		{
			res.gaps1.push_back(t_gap(0, curpos2 - curpos1));
			res.align_beg = curpos2;
		}
		else if (curpos1 > curpos2)
			res.gaps2.push_back(t_gap(0, curpos1 - curpos2));

		// number of gaps in aligned regions



		ALIGNANDUM_POSITION col_from	= src.getColFrom();
		ALIGNANDUM_POSITION col_to  	= src.getColTo();

		AlignataConstIterator it(src.begin());
		AlignataConstIterator it_end(src.end());


		ALIGNANDUM_POSITION last_col = it->mCol; ALIGNANDUM_POSITION len_col = 0;
		while (last_col < col_from)
		{
			++it;
			last_col = it->mCol;
		}
		ALIGNANDUM_POSITION last_row = it->mRow; ALIGNANDUM_POSITION len_row = 0;

		ALIGNANDUM_POSITION d_row, d_col;
		ALIGNANDUM_POSITION current_row = 0, current_col = 0;

		// start iteration at col_from + 1
		++it; len_row++; len_col++;
		for (; it != it_end; ++it)
		{

			current_row = it->mRow;
			current_col = it->mCol;

#ifdef DEBUG
//			if (verbose > LL4) cout << "current_row: " << current_row << " " << "current_col: " << current_col << endl;
#endif

			if (current_col > col_to)
				break;

			if ((d_row = current_row - last_row - 1) > 0)
			{
//				os_col << "+" << len_col << "-" << d_row;
				curpos2 += len_col;
				res.gaps2.push_back(t_gap(curpos2, d_row));

				len_col = 0;
				len_row += d_row;
			}

			if ((d_col = current_col - last_col - 1) > 0)
			{
//				os_row << "+" << len_row << "-" << d_col;
				curpos1 += len_row;
				res.gaps1.push_back(t_gap(curpos1, d_col));
				
				len_row = 0;
				len_col += d_col;
			}

			last_row = current_row;
			last_col = current_col;
			len_row++; len_col++;
		}


//		os_col << "+" << len_col << ends;
//		os_row << "+" << len_row << ends;

//
// 		last aligned region
// 
 		curpos1 += len_row;
		curpos2 += len_col;

		res.align_len = src.getLength();
		res.aligned_residues = src.getLength() - src.getNumGaps();
		
//
// 		end gap
// 
		unsigned end_seq1 =  res.len1 - curpos1;
		unsigned end_seq2 =  res.len2 - curpos2;

		if (end_seq1 > end_seq2)
			res.gaps2.push_back(t_gap(res.len2, end_seq1 - end_seq2));
		else if (end_seq1 < end_seq2)
			res.gaps1.push_back(t_gap(res.len1, end_seq2 - end_seq1));
		res.total_len = res.align_len;
		res.bitscore = src.getScore();
		
//		out_row = os_row.str(); 
//		out_col = os_col.str(); 

	}       
*/


	
}


